Préparation des données

Nous avons récolté nos données sur le site https://www.kaggle.com/datasets/sogun3/uspollution. Ces données contiennent les taux dans l’air de quatre polluants mesurés journalièrement dans plusieurs villes (stations) aux États-Unis, entre 2000 et 2016.
Nous avons restreint notre étude à 30 villes que nous avons choisi de manière arbitraire, qui correspondront ainsi à nos individus.

Importation des données

data_pollution <- read.csv("data_pollution.csv", sep=";")

data_pollution$Date.Local <- as.Date(data_pollution$Date.Local)

data_pollution <- data_pollution %>% group_by(City, Date.Local) %>% slice(1)

Dans la suite, nous avons choisi d’étudier le NO2, car ce polluant est plus pertinent dans le contexte de la qualité de l’air, en particulier dans les zones urbaines. Les émissions de NO2 sont souvent liées aux activités de combustion, telles que le trafic automobile et les installations industrielles.
L’objectif est de déterminer la tendance du taux de NO2 dans l’air aux États-Unis entre 2000 et 2016.

Description et représentations informatives des données

Statistiques descriptives

summary(data_pollution$NO2.AQI)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00   20.00   29.00   30.31   39.00  128.00
ggplot(data_pollution, aes(x = NO2.AQI)) +
  geom_histogram(binwidth = 5, fill = "#115c30", color = "black") +
  labs(title = "Distribution de NO2.AQI", x = "NO2.AQI", y = "Fréquence") +
  theme_minimal()

On observe que les valeurs de NO2 sont majoritairement comprises entre 0 et 50 ppb (parts per billion).

Tendance temporelle

ggplot(data_pollution, aes(x = Date.Local, y = NO2.AQI, group = City, color = City)) +
  geom_line() +
  labs(title = "Tendance temporelle de NO2.AQI", x = "Date", y = "NO2.AQI") +
  theme_minimal()

Cette représentation est peu lisible du fait du nombre important de données et d’individus. Pour mieux les appréhender, il serait judicieux d’appliquer un lissage à notre jeu de données, ce que nous ferons par la suite.

Pour une meilleure visualisation, nous prenons les données de la ville de Phoenix car c’est la première ville du jeu de données.

data_pho <- subset(data_pollution, City == "Phoenix")

ggplot(data_pho, aes(x = Date.Local, y = NO2.AQI)) +
  geom_line(color = "#09331a") +
  labs(title = "Tendance temporelle de NO2.AQI à Phoenix", x = "Date", y = "NO2.AQI") +
  theme_minimal()

La tendance du taux de NO2 dans l’air à Phoenix semble être à la baisse entre 2000 et 2016. Nous allons voir s’il en est de même pour toutes nos villes.

Identification des pics de pollution

#high_NO2_days <- data_pollution %>%
#  filter(NO2.AQI > 100)  # Remplacez 'seuil' par votre valeur seuil

#ggplot(data_pollution, aes(x = Date.Local, y = NO2.AQI)) +
#  geom_line(col="#09331a") +
#  geom_point(data = high_NO2_days, aes(x = Date.Local, y = NO2.AQI), col="#58a878") +
#  labs(title = "Identification des pics de pollution de NO2.AQI", x = "Date", y = "NO2.AQI") +
#  theme_minimal()

Avec le temps, on observe de moins en moins de pics de pollution, ce qui suggère une diminution globale de la pollution des villes aux États-Unis.

Représentation mensuelle

data_pollution$Month <- format(data_pollution$Date.Local, "%m")

monthly_avg <- data_pollution %>%
  group_by(Month) %>%
  summarize(Avg_NO2.AQI = mean(NO2.AQI, na.rm = TRUE))

ggplot(monthly_avg, aes(x = Month, y = Avg_NO2.AQI, group = 1)) +
  geom_line(col="#09331a") +
  geom_point(col="#09331a") +
  labs(title = "Moyenne mensuelle de NO2.AQI",
       x = "Mois",
       y = "Moyenne de NO2.AQI") +
  theme_minimal()

On observe un taux de pollution plus important l’hiver, probablement en raison d’une demande de chauffage qui augmente. De ce fait, de nombreuses personnes utilisent des combustibles fossiles tels que le charbon, le bois, le gaz naturel ou le fioul pour chauffer leurs maisons. La combustion de ces combustibles libère du NO2 dans l’atmosphère.

Lissage des données

Nous allons à présent lisser nos données à l’aide d’une méthode appropriée afin de mettre en évidence la tendance du taux de NO2 dans l’air à long terme.

Lissage avec splines

Les splines permettent de capturer des tendances non linéaires, ce qui est souvent le cas dans les séries temporelles de données environnementales comme la pollution atmosphérique. La concentration de polluants peut varier de manière non linéaire au fil du temps en raison de divers facteurs saisonniers, météorologiques, ou d’événements spécifiques.

y <- data_pho$NO2.AQI
splbasis = create.bspline.basis(rangeval=c(1, nrow(data_pho)), norder=4, breaks=seq(1, nrow(data_pho), length=100))
Phi = getbasismatrix(1:nrow(data_pho), splbasis)
chat = solve(crossprod(Phi), crossprod(Phi, y))
fhat = fd(chat, splbasis)
plot(y, type="l", cex=0.5, col="#09331a")
lines(fhat, col="#58a878", lwd=2)

chat2 = Data2fd(1:nrow(data_pho), y, basisobj = splbasis)
fhat2 = eval.fd(1:nrow(data_pho), chat2)
chatsb = smooth.basis(argvals=1:nrow(data_pho), y, fdParobj = splbasis)

Ici, pour les données de la ville de Phoenix, on observe que le lissage a permis de capter des variations saisonnière, probablement liées aux différents épisodes hivernaux qui traduisent une augmentation du taux de NO2 dans l’air en raison de la combustion liée au chauffage urbain.

Lissage avec splines pénalisées

Pour une meilleure adaptibilité à nos données, nous allons à présent appliquer un lissage avec les splines pénalisées.
Nous testons d’abord avec la ville de Phoenix pour tenter d’obtenir le meilleur lissage.

data_pho$date_num <- as.numeric(data_pho$Date.Local - min(data_pho$Date.Local))

lambda_values <- seq(0.1, 0.9, by = 0.1)

all_smoothed <- data.frame()

for (lambda in lambda_values) {
  fit <- smooth.spline(data_pho$date_num, data_pho$NO2.AQI, spar = 1-lambda)
  
  smoothed <- data.frame(
    date = data_pho$Date.Local,
    smoothed_NO2_AQI = predict(fit, data_pho$date_num)$y,
    spar = as.factor(1-lambda)
  )
  
  all_smoothed <- rbind(all_smoothed, smoothed)
}

ggplot(all_smoothed, aes(x = date, y = smoothed_NO2_AQI, color = spar)) +
  geom_line() +
  labs(title = "Comparaison des courbes lissées avec splines pénalisées",
       x = "Date",
       y = "NO2 AQI Lissé") +
  theme_minimal() +
  theme(legend.title = element_text(size = 10), legend.text = element_text(size = 8))

On n’observe pas de différence significative entre les pénalités lorsqu’elles sont inférieures à 0.5.

On aplique le lissage à toutes les villes.

data_pollution$date_num <- as.numeric(data_pollution$Date.Local - min(data_pollution$Date.Local))

lambda_value <- 0.5

all_smoothed <- data.frame()

for (city in unique(data_pollution$City)) {
  data_city <- subset(data_pollution, City == city)
  fit <- smooth.spline(data_city$date_num, data_city$NO2.AQI, spar = 1 - lambda_value)
  
  smoothed <- data.frame(
    city = city,
    date = data_city$Date.Local,
    smoothed_NO2_AQI = predict(fit, data_city$date_num)$y
  )
  
  all_smoothed <- rbind(all_smoothed, smoothed)
}

ggplot(all_smoothed, aes(x = date, y = smoothed_NO2_AQI, color = city)) +
  geom_line() +
  labs(title = "Courbes lissées avec splines pénalisées (lambda = 0.5)",
       x = "Date",
       y = "NO2 AQI Lissé") +
  theme_minimal() +
  theme(legend.title = element_text(size = 10), legend.text = element_text(size = 8))

Ce lissage nous permet donc d’observer la tendance globale du taux de NO2 dans l’air qui semble diminuer au fil du temps. Par ailleurs, on perçoit davantage les variations saisonnières sur ce dernier graphique.

Analyse fonctionnelle

L’abalyse temporelle peut nous aider à identifier les tendances temporelles significatives dans nos données lissées.

data_long <- melt(data_pollution, id.vars=c("City", "Date.Local"), variable.name="Pollutant")
data_long$Date.Local <- as.Date(data_long$Date.Local)
data_long$Pollutant <- as.factor(data_long$Pollutant)

no2_data <- data_long[data_long$Pollutant == "NO2.AQI", ]
no2_fd <- smooth.spline(no2_data$Date.Local, no2_data$value, df = 10)

plot(no2_fd, col="#09331a")

summary(no2_fd)
##            Length Class             Mode   
## x          5996   -none-            numeric
## y          5996   -none-            numeric
## w          5996   -none-            numeric
## yin        5996   -none-            numeric
## tol           1   -none-            numeric
## data          3   -none-            list   
## no.weights    1   -none-            logical
## n             1   -none-            numeric
## lev        5996   -none-            numeric
## cv            1   -none-            logical
## cv.crit       1   -none-            numeric
## pen.crit      1   -none-            numeric
## crit          1   -none-            numeric
## df            1   -none-            numeric
## spar          1   -none-            numeric
## ratio         1   -none-            numeric
## lambda        1   -none-            numeric
## iparms        5   -none-            numeric
## auxM          0   -none-            NULL   
## fit           5   smooth.spline.fit list   
## call          4   -none-            call
no2_basis <- create.bspline.basis(rangeval = range(data_pollution$Date.Local), nbasis = 10)

no2_fd <- smooth.basis(argvals = data_pollution$Date.Local, y = data_pollution$NO2.AQI, fdParobj = no2_basis)

plot(no2_fd, main = "Analyse de données fonctionnelles de NO2.AQI", col="#09331a")

## [1] "done"

La courbe obtenue à partir de l’analyse fonctionnelle peut montrer la tendance générale de la pollution de NO2.AQI au fil des années. On observe qu’il y a une diminution dans les niveaux de pollution.

Moyenne et variance fontionnelles

mean_smoothed <- all_smoothed %>%
  group_by(date) %>%
  summarise(mean_smoothed_variable = mean(smoothed_NO2_AQI),
            sd_smoothed_variable = sd(smoothed_NO2_AQI))

ggplot() +
  geom_ribbon(data = mean_smoothed, aes(x = date, ymin = mean_smoothed_variable - sd_smoothed_variable, ymax = mean_smoothed_variable + sd_smoothed_variable), fill = "#09331a", alpha = 0.3) +
  geom_line(data = mean_smoothed, aes(x = date, y = mean_smoothed_variable), color = "#09331a", size = 1) +
  labs(title = "Moyenne et Variance fonctionnelle",
       x = "Date", y = "Variable Lissée") +
  theme_minimal() +
  theme(legend.position = "bottom")
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

La ligne verte représente la moyenne fonctionnelle lissée au fil du temps. Cette ligne présente une tendance significative, qui suggère une augmentation générale de la qualité de l’air (ou diminution du taux de NO2 dans l’air) sur la période étudiée.
La bande verte délimitée par les régions ombrées représente l’intervalle de confiance autour de la moyenne lissée. Cet intervalle suggère la variabilité des données lissées à différentes périodes. Une bande plus large indique une plus grande variabilité.
La bande d’intervalle de confiance reste relativement constante sur la période, ce qui indique une stabilité de la tendance.

Préparation des données de 2015 sans valeurs manquantes pour l’ACP

Pour réaliser une ACP fonctionnelle, nous allons nous baser sur les données de 2015. Étant donné que les données présentent une certaine stabilité temporelle au fil des années, en se concentrant sur une année spécifique (par exemple, 2015), on peut simplifier l’analyse en supposant que les tendances temporelles restent relativement constantes pendant cette période.
Comme notre jeu de données est incomplet (toutes les villes/stations n’ont pas enregistré de données pour toutes les dates), nous avons choisi d’imputer les valeurs manquantes par la moyenne. C’est une méthode d’imputation simple, rapide et adaptée à un jeu de données volumineux car elle n’est pas trop coûteuses en termes de ressources computationnelles.
Malgré une tendance à la baisse du taux de NO2 dans l’air au fil du temps, nous observons une stabilité générale, ce qui peut être considéré comme un facteur atténuant. Cette stabilité temporelle contribue à minimiser le risque de biais important, renforçant ainsi la fiabilité de nos analyses.

all_smoothed$Day <- format(all_smoothed$date, "%b%d")
all_smoothed$Year <- format(all_smoothed$date, "%Y")

new_df <- all_smoothed %>%
  pivot_wider(names_from = city, values_from = smoothed_NO2_AQI)

data_2015 <- new_df[which(new_df$Year == "2015"),]
data_2015 <- column_to_rownames(data_2015, var = "Day")
data_2015 <- data_2015[, -c(1:2)]

data_2015 <- data_2015 %>% mutate_all(as.numeric)

data_2015 <- data_2015[, colSums(is.na(data_2015)) != nrow(data_2015)]

#remplacement des valeurs vides par la moyenne de la colonne
imputer <- function(col) {
  non_na_indices <- which(!is.na(col))
  for (i in which(is.na(col))) {
    col[i] <- mean(col[non_na_indices])
  }

  return(col)
}

data_2015 <- as.data.frame(apply(data_2015, 2, imputer))

data_2015 <- as.matrix(data_2015)
villes <- c("Boston", "Camden", "Charlotte", "Cleveland", "Dallas", "Houston", 
            "Kansas City", "Long Beach", "Los Angeles", "Louisville", "New York",
            "Oakland", "Oklahoma City", "Phoenix", "Portland", "Salt Lake City", 
            "Seattle", "Tucson", "Victorville", "Washington")

directions <- c("Nord-Est", "Nord-Est", "Sud-Est", "Nord-Est", "Sud",
                "Sud-Ouest", "Centre", "Sud-Ouest", "Sud-Ouest", "Centre",
                "Nord-Est", "Sud-Ouest", "Centre", "Sud-Ouest", "Nord-Ouest",
                "Centre", "Nord-Ouest", "Sud-Ouest", "Sud-Ouest", "Nord-Est")

df_direction <- data.frame(ville = villes, direction = directions)

ACP fonctionnelle

splbasis = create.bspline.basis(c(1,365),norder=4,breaks=seq(1,365,30))
fdparData = fdPar(splbasis,Lfdobj = 2,lambda=lambda)
Datasmooth = smooth.basis(1:365,data_2015,fdParobj = fdparData)

DataACPF = pca.fd(Datasmooth$fd,nharm=4,centerfns = TRUE)
DataACPF$varprop
## [1] 0.82844388 0.08714888 0.03805162 0.02183855
cumsum(DataACPF$varprop)
## [1] 0.8284439 0.9155928 0.9536444 0.9754829
plot(DataACPF$harmonics)

## [1] "done"

Ou, pour une interprétation plus aisée, la moyenne augmentée ou diminuée des premières fonctions propres.

On peut utiliser la fonction plot.pca.fd pour cette représentation.

plot.pca.fd(DataACPF)

  • Représentation des individus (les villes) dans le premier plan factoriel.

Les scores des individus sont stockés dans la matrice scores de DataACPF

head(DataACPF$scores)
##             [,1]       [,2]       [,3]       [,4]
## [1,]    7.296737   2.380958 -8.5722967 -25.653871
## [2,]   20.642756 -13.721377  5.6655690 -28.449646
## [3,] -129.498997 -26.423794 -1.4913114   2.242192
## [4,]    1.980791  53.041172  5.9861753 -37.426080
## [5,]  -96.791197 -30.270904 -0.7867414   5.230908
## [6,]  -16.025385  -6.858711 -5.2462953   9.810082
plot(DataACPF$scores[,1],DataACPF$scores[,2],pch=20,xlab="FPC 1", ylab="FPC 2",type="n")
nomsvilles = colnames(data_2015)
position <- as.factor(df_direction$direction)
text(DataACPF$scores[,1],DataACPF$scores[,2],labels=nomsvilles,cex=0.7, col = as.numeric(position))
legend("topright", legend = levels(position), col = 1:length(levels(position)), pch = 1, cex = 0.8)

On observe la formation de clusters, notamment en fonction de la position géographique sur la carte des Etats-Unis. Et que la quantité de NO2 dans l’air n’est pas uniquement liée à la saison comme on l’a indiqué précédemment, mais qu’il s’agit également d’un phénomène géographique.